Analysis of GLDS-38 from NASA GeneLab
This R markdown file was auto-generated by the iDEP website Using iDEP 0.91, originally by Steven Xijin.Ge@sdstate.edu
Ge SX, Son EW, Yao R: iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinformatics 2018, 19(1):534. PMID:30567491
First we set up the working directory to where the files are saved.
setwd('~/Documents/iDEP/iDEP-GLDS/GLDS38')
input_biclustMethod <- 'BCCC()'
R packages used
library(RSQLite,verbose=FALSE) # for database connection
library(gplots,verbose=FALSE) # for hierarchical clustering
library(ggplot2,verbose=FALSE) # graphics
library(e1071,verbose=FALSE) # computing kurtosis
library(DT,verbose=FALSE) # for renderDataTable
library(plotly,verbose=FALSE) # for interactive heatmap
library(reshape2,verbose=FALSE) # for melt correlation matrix in heatmap
# From Data Read Function
library(edgeR,verbose=FALSE) # count data D.E.
library(DESeq2,verbose=FALSE) # count data analysis, DEG.DESeq2
# TSNE Plot, tSNEgenePlot
library(Rtsne,verbose=FALSE)
# PGSA Pathway PGSEA Pathway, PGSEAplot
library(PGSEA,verbose=FALSE)
## Warning: Package 'KEGG.db' is deprecated and will be removed from Bioconductor
## version 3.12
# DEG.limma
library(limma,verbose=FALSE) # Differential expression
library(statmod,verbose=FALSE)
# enrichment plot
library(dendextend) # customizing tree
# enrich.net2, moduleNetwork
library(igraph)
# Stringdb_geneList, StringDB_GO_enrichmentData, stringDB_network1
# StringDB_network_link
library(STRINGdb,verbose=FALSE)
# gagePathwayData
library(gage,verbose=FALSE) # pathway analysis
#fgseaPathwayData
library(fgsea,verbose=FALSE) # fast GSEA
#ReactomePAPathwayData
library(ReactomePA,verbose=FALSE) # pathway analysis
# KeggImage
library(pathview)
# genomePlot, genomePlotDataPre
library(PREDA,verbose=FALSE) # showing expression on genome
library(PREDAsampledata,verbose=FALSE)
library(hgu133plus2.db,verbose=FALSE)
# biclustering
library(biclust,verbose=FALSE)
if(input_biclustMethod == "BCQU()" )
library(QUBIC,verbose=FALSE) # have trouble installing on Linux
if(input_biclustMethod == "BCUnibic()" )
library(runibic,verbose=FALSE) # Test biclustMethod dependant qubic runibic
# wgcna
library(WGCNA)
library(flashClust,verbose=FALSE)
if(file.exists('iDEP_core_functions.R'))
source('iDEP_core_functions.R') else
source('https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/iDEP_core_functions.R')
#Each row of this matrix represents a color scheme;
mycolors = sort(rainbow(20))[c(1,20,10,11,2,19,3,12,4,13,5,14,6,15,7,16,8,17,9,18)]
hmcols <- colorRampPalette(rev(c("#D73027", "#FC8D59", "#FEE090", "#FFFFBF",
"#E0F3F8", "#91BFDB", "#4575B4")))(75)
heatColors = rbind( greenred(75), bluered(75), colorpanel(75,"green","black","magenta"),colorpanel(75,"blue","yellow","red"),hmcols )
rownames(heatColors) = c("Green-Black-Red","Blue-White-Red","Green-Black-Magenta","Blue-Yellow-Red","Blue-white-brown")
We are using the downloaded gene expression file where gene IDs has been converted to Ensembl gene IDs. This is because the ID conversion database is too large to download. You can use your original file if your file uses Ensembl ID, or you do not want to use the pathway files available in iDEP (or it is not available).
inputFile_ <- 'GLDS38_Expression.csv'
sampleInfoFile_ <- 'GLDS38_Sampleinfo.csv'
gldsMetadataFile_ <- 'GLDS38_Metadata.csv'
geneInfoFile_ <- 'Arabidopsis_thaliana__athaliana_eg_gene_GeneInfo.csv' #Gene symbols, location etc.
geneSetFile_ <- 'Arabidopsis_thaliana__athaliana_eg_gene.db' # pathway database in SQL; can be GMT format
STRING10_speciesFile_ <- 'https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv'
Parameters for reading data
input_missingValue_ <- 'geneMedian' #Missing values imputation method
input_dataFileFormat_ <- 1 #1- read counts, 2 FKPM/RPKM or DNA microarray (Brake up function?)
input_minCounts_ <- 0.5 #Min counts
input_NminSamples_ <- 1 #Minimum number of samples
input_countsLogStart_ <- 4 #Pseudo count for log CPM
input_CountsTransform_ <- 1 #Methods for data transformation of counts. 1-EdgeR's logCPM 2-VST, 3-rlog
readMetadata.out_ <- readMetadata(gldsMetadataFile_)
library(knitr) # install if needed. for showing tables with kable
library(kableExtra)
kable( readMetadata.out_ ) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| FLT_Rep1 | FLT_Rep2 | FLT_Rep3 | GC_Rep1 | GC_Rep2 | GC_Rep3 | LN2_Rep1 | LN2_Rep2 | LN2_Rep3 | RNAlat_Rep1 | RNAlat_Rep2 | RNAlat_Rep3 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Sample.LongId | Atha.WT.Col.0.sl.FLT.Rep1.G1S1.RNAseq.RNAseq | Atha.WT.Col.0.sl.FLT.Rep2.G1S2.RNAseq.RNAseq | Atha.WT.Col.0.sl.FLT.Rep3.G1S3.RNAseq.RNAseq | Atha.WT.Col.0.sl.GC.Rep1.G2S1.RNAseq.RNAseq | Atha.WT.Col.0.sl.GC.Rep2.G2S2.RNAseq.RNAseq | Atha.WT.Col.0.sl.GC.Rep3.G2S3.RNAseq.RNAseq | Atha.WT.Col.0.sl.LN2.Rep1.n2.1.RNAseq.RNAseq | Atha.WT.Col.0.sl.LN2.Rep2.n2.2.RNAseq.RNAseq | Atha.WT.Col.0.sl.LN2.Rep3.n2.3.RNAseq.RNAseq | Atha.WT.Col.0.sl.RNAlat.Rep1.rl1.RNAseq.RNAseq | Atha.WT.Col.0.sl.RNAlat.Rep2.rl2.RNAseq.RNAseq | Atha.WT.Col.0.sl.RNAlat.Rep3.rl3.RNAseq.RNAseq |
| Sample.Id | Atha.WT.Col.0.sl.FLT.Rep1.G1S1 | Atha.WT.Col.0.sl.FLT.Rep2.G1S2 | Atha.WT.Col.0.sl.FLT.Rep3.G1S3 | Atha.WT.Col.0.sl.GC.Rep1.G2S1 | Atha.WT.Col.0.sl.GC.Rep2.G2S2 | Atha.WT.Col.0.sl.GC.Rep3.G2S3 | Atha.WT.Col.0.sl.LN2.Rep1.n2.1 | Atha.WT.Col.0.sl.LN2.Rep2.n2.2 | Atha.WT.Col.0.sl.LN2.Rep3.n2.3 | Atha.WT.Col.0.sl.RNAlat.Rep1.rl1 | Atha.WT.Col.0.sl.RNAlat.Rep2.rl2 | Atha.WT.Col.0.sl.RNAlat.Rep3.rl3 |
| Sample.Name | Atha_WT-Col-0_sl_FLT_Rep1_G1S1 | Atha_WT-Col-0_sl_FLT_Rep2_G1S2 | Atha_WT-Col-0_sl_FLT_Rep3_G1S3 | Atha_WT-Col-0_sl_GC_Rep1_G2S1 | Atha_WT-Col-0_sl_GC_Rep2_G2S2 | Atha_WT-Col-0_sl_GC_Rep3_G2S3 | Atha_WT-Col-0_sl_LN2_Rep1_n2-1 | Atha_WT-Col-0_sl_LN2_Rep2_n2-2 | Atha_WT-Col-0_sl_LN2_Rep3_n2-3 | Atha_WT-Col-0_sl_RNAlat_Rep1_rl1 | Atha_WT-Col-0_sl_RNAlat_Rep2_rl2 | Atha_WT-Col-0_sl_RNAlat_Rep3_rl3 |
| GLDS | 38 | 38 | 38 | 38 | 38 | 38 | 38 | 38 | 38 | 38 | 38 | 38 |
| Accession | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 | GLDS-38 |
| Hardware | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC | BRIC |
| Tissue | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling | Etiolated seedling |
| Age | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days | 8 days |
| Organism | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana | Arabidopsis thaliana |
| Ecotype | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 | Col-0 |
| Genotype | WT | WT | WT | WT | WT | WT | WT | WT | WT | WT | WT | WT |
| Variety | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT | Col-0 WT |
| Radiation | Cosmic radiation | Cosmic radiation | Cosmic radiation | Background Earth | Background Earth | Background Earth | Background Earth | Background Earth | Background Earth | Background Earth | Background Earth | Background Earth |
| Gravity | Microgravity | Microgravity | Microgravity | Terrestrial | Terrestrial | Terrestrial | Terrestrial | Terrestrial | Terrestrial | Terrestrial | Terrestrial | Terrestrial |
| Developmental | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings | Etiolated 8 day old seedlings |
| Time.series.or.Concentration.gradient | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point | Single time point |
| Light | Dark | Dark | Dark | Dark | Dark | Dark | Dark | Dark | Dark | Dark | Dark | Dark |
| Assay..RNAseq. | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling | RNAseq Transcription and Proteomic Profiling |
| Temperature | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 | 22-24 |
| Treatment.type | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity | Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity |
| Treatment.intensity | X | X | X | X | X | X | X | X | X | X | X | X |
| Treament.timing | X | X | X | X | X | X | X | X | X | X | X | X |
| Preservation.Method. | RNAlater | RNAlater | RNAlater | RNAlater | RNAlater | RNAlater | Liquid Nitrogen | Liquid Nitrogen | Liquid Nitrogen | RNAlater | RNAlater | RNAlater |
readData.out <- readData(inputFile_, input_missingValue=input_missingValue_, input_dataFileFormat = input_dataFileFormat_, input_minCounts=input_minCounts_, input_NminSamples = input_NminSamples_, input_countsLogStart=input_countsLogStart_, input_CountsTransform = input_CountsTransform_)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
kable( head(readData.out$data) ) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| FLT_Rep1 | FLT_Rep2 | FLT_Rep3 | GC_Rep1 | GC_Rep2 | GC_Rep3 | LN2_Rep1 | LN2_Rep2 | LN2_Rep3 | RNAlat_Rep1 | RNAlat_Rep2 | RNAlat_Rep3 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ATCG00490 | 19.02600 | 20.20307 | 20.47902 | 19.35267 | 20.26267 | 19.81133 | 17.85944 | 18.14881 | 18.28268 | 19.90399 | 19.24373 | 18.83189 |
| AT2G41310 | 18.59846 | 19.97237 | 19.85194 | 20.18207 | 18.98753 | 19.11273 | 18.29485 | 18.20631 | 18.20418 | 18.69872 | 18.47971 | 18.18740 |
| ATCG00020 | 17.76784 | 18.79161 | 19.31738 | 17.92858 | 19.00300 | 18.51242 | 17.11173 | 17.35404 | 17.37928 | 18.97660 | 18.68141 | 18.13132 |
| AT3G21720 | 17.38043 | 17.29131 | 17.40174 | 18.45645 | 16.52706 | 17.41765 | 16.35671 | 15.58054 | 15.69393 | 14.72961 | 14.46187 | 14.90556 |
| AT2G07671 | 16.80806 | 17.55228 | 17.74476 | 16.89524 | 17.38578 | 17.32608 | 15.95185 | 15.89930 | 15.66004 | 16.40474 | 16.05359 | 15.88383 |
| ATCG00280 | 16.09754 | 17.02548 | 17.38910 | 17.07244 | 17.48614 | 17.20404 | 15.47935 | 15.34283 | 15.45231 | 17.24330 | 16.73648 | 16.46234 |
readSampleInfo.out <- readSampleInfo(sampleInfoFile_, readData.out)
kable( readSampleInfo.out ) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| Gravity | Preservation.Method. | |
|---|---|---|
| FLT_Rep1 | Microgravity | RNAlater |
| FLT_Rep2 | Microgravity | RNAlater |
| FLT_Rep3 | Microgravity | RNAlater |
| GC_Rep1 | Terrestrial | RNAlater |
| GC_Rep2 | Terrestrial | RNAlater |
| GC_Rep3 | Terrestrial | RNAlater |
| LN2_Rep1 | Terrestrial | Liquid Nitrogen |
| LN2_Rep2 | Terrestrial | Liquid Nitrogen |
| LN2_Rep3 | Terrestrial | Liquid Nitrogen |
| RNAlat_Rep1 | Terrestrial | RNAlater |
| RNAlat_Rep2 | Terrestrial | RNAlater |
| RNAlat_Rep3 | Terrestrial | RNAlater |
input_selectOrg_ ="NEW"
input_selectGO_ <- 'GOBP' #Gene set category
input_noIDConversion_ = TRUE
allGeneInfo.out_ <- geneInfo(geneInfoFile_)
converted.out_ = NULL
convertedData.out <- convertedData(converted.out=converted.out_, readData.out=readData.out, input_noIDConversion=input_noIDConversion_)
nGenesFilter(readData.out=readData.out, converted.out=converted.out_, convertedData.out=convertedData.out, input_noIDConversion=input_noIDConversion_)
## [1] "16156 genes in 12 samples. 16117 genes passed filter.\n Original gene IDs used."
convertedCounts.out_ <- convertedCounts(readData.out=readData.out, converted.out=converted.out_) # converted counts, just for compatibility
# Read counts per library
parDefault = par()
par(mar=c(12,4,2,2))
# barplot of total read counts
x <- readData.out$rawCounts
groups = as.factor( detectGroups(colnames(x ) ) )
if(nlevels(groups)<=1 | nlevels(groups) >20 )
col1 = 'green' else
col1 = rainbow(nlevels(groups))[ groups ]
barplot( colSums(x)/1e6,
col=col1,las=3, main="Total read counts (millions)")
readCountsBias(readData.out=readData.out, readSampleInfo.out=readSampleInfo.out) # detecting bias in sequencing depth
## [1] 0.02579404
## [1] 0.1547212
## [1] 0.05122098
## [1] "Warning! Sequencing depth bias detected. Total read counts are significantly different among sample groups (p= 2.58e-02 ) based on ANOVA."
# Box plot
x = readData.out$data
boxplot(x, las = 2, col=col1,
ylab='Transformed expression levels',
main='Distribution of transformed data')
#Density plot
par(parDefault)
## Warning in par(parDefault): graphical parameter "cin" cannot be set
## Warning in par(parDefault): graphical parameter "cra" cannot be set
## Warning in par(parDefault): graphical parameter "csi" cannot be set
## Warning in par(parDefault): graphical parameter "cxy" cannot be set
## Warning in par(parDefault): graphical parameter "din" cannot be set
## Warning in par(parDefault): graphical parameter "page" cannot be set
densityPlot(readData.out = readData.out, mycolors)
# Scatter plot of the first two samples
plot(x[,1:2],xlab=colnames(x)[1],ylab=colnames(x)[2],
main='Scatter plot of first two samples')
####plot gene or gene family
input_selectOrg_ ="BestMatch"
input_geneSearch_ <- 'HOXA' #Gene ID for searching
genePlot(allGeneInfo.out_, convertedData.out=convertedData.out, input_selectOrg=input_selectOrg_, input_geneSearch=input_geneSearch_)
## NULL
input_useSD_ <- 'FALSE' #Use standard deviation instead of standard error in error bar?
geneBarPlotError(allGeneInfo.out_, input_selectOrg_, input_geneSearch_, input_useSD_)
## NULL
# hierarchical clustering tree
x <- readData.out$data
maxGene <- apply(x,1,max)
# remove bottom 25% lowly expressed genes, which inflate the PPC
x <- x[which(maxGene > quantile(maxGene)[1] ) ,]
plot(as.dendrogram(hclust2( dist2(t(x)))), ylab="1 - Pearson C.C.", type = "rectangle")
#Correlation matrix
input_labelPCC_ <- TRUE #Show correlation coefficient?
correlationMatrix(input_labelPCC = input_labelPCC_)
# Parameters for heatmap
input_nGenes_ <- 1000 #Top genes for heatmap
input_geneCentering_ <- TRUE #centering genes ?
input_sampleCentering_ <- FALSE #Center by sample?
input_geneNormalize_ <- FALSE #Normalize by gene?
input_sampleNormalize_ <- FALSE #Normalize by sample?
input_noSampleClustering_ <- FALSE #Use original sample order
input_heatmapCutoff_ <- 4 #Remove outliers beyond number of SDs
input_distFunctions_ <- 1 #which distant funciton to use
input_hclustFunctions_ <- 1 #Linkage type
input_heatColors1_ <- 1 #Colors
input_selectFactorsHeatmap_ <- 'Gravity' #Sample coloring factors
png('heatmap.png', width = 10, height = 15, units = 'in', res = 300)
staticHeatmap(heatColors, input_nGenes_, input_geneCentering_, input_sampleCentering_, input_geneNormalize_, input_sampleNormalize_, input_noSampleClustering=input_noSampleClustering_, input_heatmapCutoff=input_heatmapCutoff_, input_distFunctions=input_distFunctions_, input_hclustFunctions=input_hclustFunctions_, input_heatColors1=input_heatColors1_, input_selectFactorsHeatmap=input_selectFactorsHeatmap_)
dev.off()
## png
## 2
[heatmap] (heatmap.png)
heatmapPlotly(heatColors, allGeneInfo.out=allGeneInfo.out_, input_geneCentering=input_geneCentering_, input_sampleCentering = input_sampleCentering_, input_geneNormalize = input_geneNormalize_, input_sampleNormalize = input_sampleNormalize_, input_heatColors1=input_heatColors1_) # interactive heatmap using Plotly
input_nGenesKNN_ <- 2000 #Number of genes fro k-Means
input_nClusters_ <- 4 #Number of clusters
maxGeneClustering_ = 12000
input_kmeansNormalization_ <- 'geneMean' #Normalization
input_KmeansReRun_ <- 0 #Random seed
distributionSD(input_nGenesKNN_) #Distribution of standard deviations
KmeansNclusters(input_nGenesKNN_) #Number of clusters
Kmeans.out = Kmeans(maxGeneClustering=maxGeneClustering_, input_nGenesKNN=input_nGenesKNN_, input_nClusters=input_nClusters_, input_kmeansNormalization=input_kmeansNormalization_, input_KmeansReRun=input_KmeansReRun_) #Running K-means
KmeansHeatmap(.mycolors = mycolors, .heatColors = heatColors, .input_heatColors1=input_heatColors1_) #Heatmap for k-Means
#Read gene sets for enrichment analysis
sqlite_ <- dbDriver('SQLite') #sqlite changed to default in most functions
input_selectGO3_ <- 'GOBP' #Gene set category
input_minSetSize_ <- 15 #Min gene set size
input_maxSetSize_ <- 2000 #Max gene set size
GeneSets.out <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO3_,input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
# Alternatively, users can use their own GMT files by
#GeneSets.out <- readGMTRobust('somefile.GMT')
results <- KmeansGO(input_nClusters=input_nClusters_) #Enrichment analysis for k-Means clusters
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| Cluster | adj.Pval | Genes | Pathways |
|---|---|---|---|
| A | 7.99e-122 | 219 | Organonitrogen compound biosynthetic process |
| 2.84e-118 | 155 | Translation | |
| 4.26e-118 | 155 | Peptide biosynthetic process | |
| 1.81e-116 | 159 | Amide biosynthetic process | |
| 3.07e-113 | 156 | Peptide metabolic process | |
| 1.22e-106 | 160 | Cellular amide metabolic process | |
| 5.05e-51 | 157 | Response to abiotic stimulus | |
| 4.90e-42 | 57 | Photosynthesis | |
| 5.59e-35 | 50 | Response to cytokinin | |
| 7.30e-34 | 63 | Generation of precursor metabolites and energy | |
| B | 9.12e-48 | 81 | Cell wall organization or biogenesis |
| 9.12e-48 | 73 | Cell wall organization | |
| 2.49e-47 | 74 | External encapsulating structure organization | |
| 1.20e-40 | 128 | Response to abiotic stimulus | |
| 2.58e-36 | 74 | Drug metabolic process | |
| 8.64e-36 | 111 | Small molecule metabolic process | |
| 2.30e-34 | 82 | Carbohydrate metabolic process | |
| 2.27e-32 | 77 | Response to inorganic substance | |
| 6.60e-28 | 42 | Plant-type cell wall organization or biogenesis | |
| 1.13e-27 | 51 | Polysaccharide metabolic process | |
| C | 3.98e-50 | 79 | Response to external stimulus |
| 1.34e-42 | 63 | Response to external biotic stimulus | |
| 1.34e-42 | 63 | Response to other organism | |
| 2.33e-42 | 63 | Response to biotic stimulus | |
| 4.56e-40 | 42 | Secondary metabolic process | |
| 1.25e-38 | 63 | Defense response | |
| 2.03e-38 | 67 | Multi-organism process | |
| 2.35e-35 | 78 | Response to abiotic stimulus | |
| 4.27e-34 | 26 | Indole-containing compound metabolic process | |
| 1.99e-33 | 49 | Defense response to other organism | |
| D | 1.20e-50 | 151 | Response to abiotic stimulus |
| 1.95e-47 | 127 | Response to oxygen-containing compound | |
| 2.73e-46 | 98 | Response to inorganic substance | |
| 1.10e-39 | 101 | Response to acid chemical | |
| 2.11e-39 | 113 | Cellular response to chemical stimulus | |
| 6.77e-37 | 126 | Response to organic substance | |
| 1.17e-31 | 107 | Response to hormone | |
| 4.57e-31 | 107 | Response to endogenous stimulus | |
| 3.63e-27 | 75 | Response to external biotic stimulus | |
| 3.63e-27 | 75 | Response to other organism |
input_seedTSNE_ <- 0 #Random seed for t-SNE
input_colorGenes_ <- TRUE #Color genes in t-SNE plot?
tSNEgenePlot(input_seedTSNE=input_seedTSNE_, input_colorGenes=input_colorGenes_) #Plot genes using t-SNE
input_selectFactors_ <- 'Gravity' #Factor coded by color
input_selectFactors2_ <- 'Preservation.Method.' #Factor coded by shape
input_tsneSeed2_ <- 0 #Random seed for t-SNE
#PCA, MDS and t-SNE plots
PCAplot(input_selectFactors_, input_selectFactors2_)
MDSplot(input_selectFactors_, input_selectFactors2_)
tSNEplot(input_selectFactors_, input_selectFactors2_, input_tsneSeed2_)
#Read gene sets for pathway analysis using PGSEA on principal components
input_selectGO6_ <- 'GOBP'
GeneSets.out_ <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO6_,input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
PCApathway(convertedData.out=convertedData.out, GeneSets.out=GeneSets.out_) # Run PGSEA analysis
cat( PCA2factor(readData.out, readSampleInfo.out=readSampleInfo.out) ) #The correlation between PCs with factors
##
## Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Gravity (p=3.91e-02).
input_CountsDEGMethod_ <- 3 #DESeq2= 3,limma-voom=2,limma-trend=1
input_limmaPval_ <- 0.1 #FDR cutoff
input_limmaFC_ <- 2 #Fold-change cutoff
input_selectModelComprions_ <- 'Gravity: Microgravity vs. Terrestrial' #Selected comparisons
input_selectFactorsModel_ <- 'Gravity' #Selected comparisons
input_selectInteractions_ <- NULL #Selected comparisons
input_selectBlockFactorsModel_ <- NULL #Selected comparisons
factorReferenceLevels.out_ <- c('Gravity:Terrestrial')
limma.out <- limma(input_dataFileFormat = input_dataFileFormat_, input_countsLogStart = input_countsLogStart_, convertedCounts.out=convertedCounts.out_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, input_selectInteractions=input_selectInteractions_, input_selectBlockFactorsModel=input_selectBlockFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DEG.data.out <- DEG.data(allGeneInfo.out_)
limma.out$comparisons
## [1] "Microgravity-Terrestrial"
input_selectComparisonsVenn_ = limma.out$comparisons[1:3] # use first three comparisons
input_UpDownRegulated_ <- FALSE #Split up and down regulated genes
vennPlot(input_selectComparisonsVenn_, input_UpDownRegulated_) # Venn diagram
sigGeneStats() # number of DEGs as figure
sigGeneStatsTable() # number of DEGs as table
## Comparisons Up Down
## Microgravity-Terrestrial Microgravity-Terrestrial 353 1071
input_selectContrast_ <- 'Microgravity-Terrestrial' #Selected comparisons
selectedHeatmap.data.out <- selectedHeatmap.data(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast=input_selectContrast_)
selectedHeatmap(.mycolors=mycolors,.heatColors=heatColors) # heatmap for DEGs in selected comparison
# Save gene lists and data into files
write.csv( selectedHeatmap.data(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat = input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast=input_selectContrast_)$genes, 'heatmap.data.csv')
write.csv(DEG.data(allGeneInfo.out_),'DEG.data.csv' )
write(AllGeneListsGMT() ,'AllGeneListsGMT.gmt')
input_selectGO2_ <- 'GOBP' #Gene set category
geneListData.out <- geneListData(allGeneInfo.out_, input_selectGO2_, input_selectOrg_, input_limmaPval_, input_limmaFC_, input_selectContrast_)
volcanoPlot(input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectContrast=input_selectContrast_)
scatterPlot(input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel=input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast_)
MAplot(.converted.out=converted.out_, .readData.out=readData.out, .input_noIDConversion=input_noIDConversion_, input_dataFileFormat=input_dataFileFormat_, input_CountsDEGMethod=input_CountsDEGMethod_, input_limmaPval=input_limmaPval_, input_limmaFC=input_limmaFC_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel_, factorReferenceLevels.out_, input_selectContrast_)
geneListGOTable.out <- geneListGOTable()
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO2_,input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
input_removeRedudantSets_ <- TRUE #Remove highly redundant gene sets?
results <- geneListGO(input_removeRedudantSets_) #Enrichment analysis
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| Direction | adj.Pval | nGenes | Pathways |
|---|---|---|---|
| Down regulated | 4.1e-37 | 143 | Response to external stimulus |
| 4.6e-36 | 117 | Response to biotic stimulus | |
| 7.5e-35 | 114 | Response to external biotic stimulus | |
| 7.5e-35 | 114 | Response to other organism | |
| 1.1e-34 | 124 | Defense response | |
| 2.3e-33 | 134 | Multi-organism process | |
| 1.7e-28 | 60 | Secondary metabolic process | |
| 1.1e-27 | 79 | Response to drug | |
| 1.5e-26 | 137 | Response to oxygen-containing compound | |
| 2.8e-26 | 156 | Response to organic substance | |
| Up regulated | 1.9e-18 | 29 | Photosynthesis |
| 5.7e-15 | 61 | Organonitrogen compound biosynthetic process | |
| 7.5e-15 | 41 | Peptide metabolic process | |
| 5.2e-14 | 38 | Translation | |
| 5.2e-14 | 38 | Peptide biosynthetic process | |
| 1.6e-13 | 42 | Cellular amide metabolic process | |
| 1.6e-13 | 39 | Amide biosynthetic process | |
| 3.1e-13 | 19 | Pigment biosynthetic process | |
| 1.1e-12 | 20 | Pigment metabolic process | |
| 1.3e-11 | 23 | Plastid organization |
STRING-db API access. We need to find the taxonomy id of your species, this used by STRING. First we try to guess the ID based on iDEP’s database. Users can also skip this step and assign NCBI taxonomy id directly by findTaxonomyID.out = 10090 # mouse 10090, human 9606 etc.
STRING10_species_ = read.csv(STRING10_speciesFile_)
ix = grep('Arabidopsis thaliana', STRING10_species_$official_name )
findTaxonomyID.out <- STRING10_species_[ix,1] # find taxonomyID
findTaxonomyID.out
## [1] 3702
Enrichment analysis using STRING
STRINGdb_geneList.out <- STRINGdb_geneList() #convert gene lists
## Warning: we couldn't map to STRING 0% of your identifiers
input_STRINGdbGO_ <- 'Process' #'Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro'
results <- stringDB_GO_enrichmentData(input_STRINGdbGO_) # enrichment using STRING
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| “No significant enrichment found.” | adj.Pval |
|---|---|
| No significant enrichment found. | NULL |
PPI network retrieval and analysis
input_nGenesPPI_ <- 100 #Number of top genes for PPI retrieval and analysis
stringDB_network1(1, input_nGenesPPI_) #Show PPI network
Generating interactive PPI
write(stringDB_network_link(input_nGenesPPI_), 'PPI_results.html') # write results to html file
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")
## Warning: we couldn't map to STRING 0% of your identifiers
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")
browseURL('PPI_results.html') # open in browser
input_selectContrast1_ <- 'Microgravity-Terrestrial' #select Comparison
#input_selectContrast1 = limma.out$comparisons[3] # manually set
input_selectGO_ <- 'GOBP' #Gene set category
#input_selectGO='custom' # if custom gmt file
input_minSetSize_ <- 15 #Min size for gene set
input_maxSetSize_ <- 2000 #Max size for gene set
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO_,input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
input_pathwayPvalCutoff_ <- 0.2 #FDR cutoff
input_nPathwayShow_ <- 30 #Top pathways to show
input_absoluteFold_ <- FALSE #Use absolute values of fold-change?
input_GenePvalCutoff_ <- 1 #FDR to remove genes
input_pathwayMethod_ = 1 # 1 GAGE
gagePathwayData.out <- gagePathwayData(input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff_, input_nPathwayShow_, input_absoluteFold_, input_GenePvalCutoff_) # pathway analysis using GAGE
results <- gagePathwayData.out #Enrichment analysis for k-Means clusters
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| Direction | GAGE analysis: Microgravity vs Terrestrial | statistic | Genes | adj.Pval |
|---|---|---|---|---|
| Down | Response to chitin | -7.7014 | 104 | 1.9e-10 |
| Secondary metabolic process | -7.4839 | 265 | 1.9e-10 | |
| Response to drug | -7.3469 | 481 | 1.9e-10 | |
| Response to organonitrogen compound | -6.3515 | 212 | 1.3e-07 | |
| Secondary metabolite biosynthetic process | -6.3149 | 118 | 2.4e-07 | |
| Response to nitrogen compound | -6.2461 | 271 | 1.6e-07 | |
| Up | Photosynthesis | 10.7757 | 223 | 3.7e-21 |
| Ribosome biogenesis | 10.2165 | 342 | 6.5e-20 | |
| Ribonucleoprotein complex biogenesis | 10.0934 | 437 | 6.5e-20 | |
| NcRNA metabolic process | 9.4896 | 425 | 8.0e-18 | |
| Plastid organization | 9.039 | 257 | 6.2e-16 | |
| NcRNA processing | 8.7867 | 357 | 2.8e-15 | |
| RRNA processing | 8.2838 | 239 | 2.6e-13 | |
| RRNA metabolic process | 8.2747 | 244 | 2.6e-13 | |
| Photosynthesis, light reaction | 7.9994 | 119 | 7.6e-12 | |
| Chloroplast organization | 7.4927 | 198 | 4.2e-11 | |
| Generation of precursor metabolites and energy | 7.2638 | 391 | 7.9e-11 | |
| RNA modification | 7.1316 | 321 | 2.8e-10 | |
| Thylakoid membrane organization | 6.8664 | 46 | 1.5e-07 | |
| Tetrapyrrole biosynthetic process | 6.0669 | 70 | 7.8e-07 | |
| Ribosome assembly | 6.0176 | 76 | 1.3e-06 | |
| Ribosomal large subunit biogenesis | 5.9801 | 99 | 1.0e-06 | |
| Tetrapyrrole metabolic process | 5.964 | 93 | 7.8e-07 | |
| Porphyrin-containing compound biosynthetic process | 5.9033 | 67 | 1.3e-06 | |
| Porphyrin-containing compound metabolic process | 5.8821 | 92 | 1.1e-06 | |
| Protein transmembrane transport | 5.7744 | 112 | 1.4e-06 | |
| Ribonucleoprotein complex subunit organization | 5.7385 | 181 | 1.3e-06 | |
| Nucleoside monophosphate metabolic process | 5.6492 | 193 | 1.6e-06 | |
| Ribonucleoprotein complex assembly | 5.6323 | 173 | 1.8e-06 | |
| TRNA metabolic process | 5.6148 | 170 | 1.8e-06 |
pathwayListData.out_ = pathwayListData(allGeneInfo.out_, input_selectOrg_, input_selectGO_, input_pathwayMethod_)
enrichmentPlot(pathwayListData.out_, 25 )
enrichmentNetwork(pathwayListData.out_ )
enrichmentNetworkPlotly(pathwayListData.out_)
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
input_pathwayMethod_ = 3 # 1 fgsea
fgseaPathwayData.out <- fgseaPathwayData(input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff=input_pathwayPvalCutoff_, input_absoluteFold=input_absoluteFold_, input_nPathwayShow=input_nPathwayShow_, input_GenePvalCutoff=input_GenePvalCutoff_) #Pathway analysis using fgsea
## Warning in fgsea(pathways = gmt, stats = fold, minSize = input_minSetSize, :
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in fgseaSimple(...): There were 3 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
results <- fgseaPathwayData.out #Enrichment analysis for k-Means clusters
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| Direction | GSEA analysis: Microgravity vs Terrestrial | NES | Genes | adj.Pval |
|---|---|---|---|---|
| Up | Photosynthesis | 3.2975 | 223 | 1.2e-02 |
| Photosynthesis, light reaction | 3.0161 | 119 | 6.3e-03 | |
| Thylakoid membrane organization | 2.9239 | 46 | 4.1e-03 | |
| Translation | 2.9137 | 619 | 1.1e-01 | |
| Peptide biosynthetic process | 2.9053 | 622 | 1.1e-01 | |
| Plastid organization | 2.898 | 257 | 1.4e-02 | |
| Amide biosynthetic process | 2.8546 | 683 | 1.8e-01 | |
| Ribosome biogenesis | 2.8254 | 342 | 2.3e-02 | |
| Photosynthetic electron transport chain | 2.8113 | 46 | 4.1e-03 | |
| Chloroplast organization | 2.7408 | 198 | 9.8e-03 | |
| Ribosome assembly | 2.7189 | 76 | 4.8e-03 | |
| Tetrapyrrole biosynthetic process | 2.6943 | 70 | 4.7e-03 | |
| RRNA processing | 2.6919 | 239 | 1.3e-02 | |
| Plastid membrane organization | 2.6845 | 49 | 4.2e-03 | |
| RRNA metabolic process | 2.6757 | 244 | 1.3e-02 | |
| Chlorophyll biosynthetic process | 2.6734 | 58 | 4.3e-03 | |
| Protein localization to chloroplast | 2.6602 | 45 | 4.1e-03 | |
| Porphyrin-containing compound biosynthetic process | 2.6383 | 67 | 4.6e-03 | |
| Ribonucleoprotein complex biogenesis | 2.6344 | 437 | 3.9e-02 | |
| Ribosomal large subunit biogenesis | 2.6273 | 99 | 5.6e-03 | |
| Tetrapyrrole metabolic process | 2.6198 | 93 | 5.4e-03 | |
| Protein targeting to chloroplast | 2.6126 | 43 | 4.1e-03 | |
| Establishment of protein localization to chloroplast | 2.6126 | 43 | 4.1e-03 | |
| Chlorophyll metabolic process | 2.6057 | 81 | 4.9e-03 | |
| Porphyrin-containing compound metabolic process | 2.5984 | 92 | 5.4e-03 | |
| Peptide metabolic process | 2.5698 | 679 | 1.7e-01 | |
| Chloroplast RNA processing | 2.5648 | 19 | 3.9e-03 | |
| Pigment biosynthetic process | 2.5423 | 129 | 6.8e-03 | |
| NcRNA processing | 2.499 | 357 | 2.4e-02 | |
| Chloroplast rRNA processing | 2.4985 | 18 | 3.9e-03 |
pathwayListData.out_ = pathwayListData(allGeneInfo.out_, input_selectOrg=input_selectOrg_, input_selectGO=input_selectGO_, input_pathwayMethod_)
enrichmentPlot(pathwayListData.out_, 25, mycolors)
enrichmentNetwork(pathwayListData.out_ )
enrichmentNetworkPlotly(pathwayListData.out_)
PGSEAplot(input_selectOrg = input_selectOrg_, input_dataFileFormat=input_dataFileFormat_, input_selectGO=input_selectGO_, input_minSetSize=input_minSetSize_, input_maxSetSize=input_maxSetSize_, input_CountsDEGMethod=input_CountsDEGMethod_, input_selectModelComprions=input_selectModelComprions_, input_selectFactorsModel_, factorReferenceLevels.out=factorReferenceLevels.out_, input_selectContrast1=input_selectContrast1_, input_pathwayPvalCutoff=input_pathwayPvalCutoff_, input_nPathwayShow_) # pathway analysis using PGSEA
##
## Computing P values using ANOVA
input_selectContrast2_ <- 'Microgravity-Terrestrial' #select Comparison
#input_selectContrast2 = limma.out$comparisons[3] # manually set
input_limmaPvalViz_ <- 0.1 #FDR to filter genes
input_limmaFCViz_ <- 2 #FDR to filter genes
genomePlotly(allGeneInfo.out_, input_selectContrast2_, input_limmaPvalViz_, input_limmaFCViz_) # shows fold-changes on the genome
## Warning in eval(quote(list(...)), env): NAs introduced by coercion
## Warning in genomePlotly(allGeneInfo.out_, input_selectContrast2_,
## input_limmaPvalViz_, : NAs introduced by coercion
input_nGenesBiclust_ <- 1000 #Top genes for biclustering
input_biclustMethod_ <- 'BCCC()' #Method: 'BCCC', 'QUBIC', 'runibic' ...
biclustering.out = biclustering(input_nGenesBiclust_, input_biclustMethod_ ) # run analysis
input_selectBicluster_ <- 1 #select a cluster
biclustHeatmap(heatColors, input_heatColors1=input_heatColors1_, input_selectBicluster=input_selectBicluster_) # heatmap for selected cluster
input_selectGO4_ <- 'GOBP' #Gene set category
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO4_, input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
results <- geneListBclustGO(2, input_selectBicluster=input_selectBicluster_) #Enrichment analysis for k-Means clusters
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| adj.Pval | Genes | Pathways |
|---|---|---|
| 2.4e-134 | 308 | Response to abiotic stimulus |
| 5.9e-86 | 170 | Response to inorganic substance |
| 1.1e-66 | 216 | Response to organic substance |
| 5.7e-64 | 182 | Oxidation-reduction process |
| 1.3e-63 | 197 | Organonitrogen compound biosynthetic process |
| 2.2e-63 | 173 | Response to external stimulus |
| 4.5e-62 | 144 | Response to external biotic stimulus |
| 4.5e-62 | 144 | Response to other organism |
| 2.6e-61 | 144 | Response to biotic stimulus |
| 3.1e-61 | 105 | Response to metal ion |
input_mySoftPower_ <- 5 #SoftPower to cutoff
input_nGenesNetwork_ <- 1000 #Number of top genes
input_minModuleSize_ <- 20 #Module size minimum
wgcna.out = wgcna(2000, input_mySoftPower=input_mySoftPower_, input_nGenesNetwork=input_nGenesNetwork_, input_minModuleSize=input_minModuleSize_) # run WGCNA
## Warning: executing %dopar% sequentially: no parallel backend registered
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.6880 3.020 0.962 418.00 426.00 557.0
## 2 2 0.5030 1.310 0.949 240.00 241.00 388.0
## 3 3 0.2260 0.549 0.896 159.00 155.00 294.0
## 4 4 0.0337 0.163 0.868 113.00 108.00 233.0
## 5 5 0.0639 -0.187 0.865 85.00 79.00 191.0
## 6 6 0.2370 -0.350 0.875 66.30 60.10 159.0
## 7 7 0.4190 -0.473 0.908 53.20 46.40 135.0
## 8 8 0.5660 -0.576 0.946 43.70 37.00 117.0
## 9 9 0.6550 -0.661 0.944 36.50 30.10 103.0
## 10 10 0.7000 -0.706 0.948 30.90 24.80 91.1
## 11 12 0.7590 -0.822 0.928 23.00 17.50 73.4
## 12 14 0.7970 -0.918 0.922 17.70 12.80 60.7
## 13 16 0.8360 -0.968 0.932 14.00 9.68 51.3
## 14 18 0.8470 -1.020 0.929 11.40 7.55 44.1
## 15 20 0.8340 -1.050 0.905 9.39 5.96 38.5
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
softPower(wgcna.out) # soft power curve
modulePlot(wgcna.out) # plot modules
listWGCNA.Modules.out = listWGCNA.Modules() #modules
input_selectGO5_ <- 'GOBP' #Gene set category
# Read pathway data again
GeneSets.out <-readGeneSets( geneSetFile_,
convertedData.out, input_selectGO5_,input_selectOrg_,
c(input_minSetSize_, input_maxSetSize_) )
input_selectWGCNA.Module_ <- 'Entire network' #Select a module
input_topGenesNetwork_ <- 10 #SoftPower to cutoff
input_edgeThreshold_ <- 0.4 #Number of top genes
moduleNetwork(input_noIDConversion=input_noIDConversion_, allGeneInfo.out=allGeneInfo.out_, input_selectOrg_, input_selectWGCNA.Module_, input_topGenesNetwork_, input_edgeThreshold_, input_selectGO5_) # show network of top genes in selected module
## softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
## ..calculating connectivities..
input_removeRedudantSets_ <- TRUE #Remove redundant gene sets
results <- networkModuleGO(2, input_selectWGCNA.Module_) #Enrichment analysis of selected module
results$adj.Pval <- format( results$adj.Pval,digits=3 )
kable( results, row.names=FALSE) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%")
| adj.Pval | Genes | Pathways |
|---|---|---|
| 2.4e-134 | 308 | Response to abiotic stimulus |
| 5.9e-86 | 170 | Response to inorganic substance |
| 1.1e-66 | 216 | Response to organic substance |
| 5.7e-64 | 182 | Oxidation-reduction process |
| 1.3e-63 | 197 | Organonitrogen compound biosynthetic process |
| 2.2e-63 | 173 | Response to external stimulus |
| 4.5e-62 | 144 | Response to external biotic stimulus |
| 4.5e-62 | 144 | Response to other organism |
| 2.6e-61 | 144 | Response to biotic stimulus |
| 3.1e-61 | 105 | Response to metal ion |